package trim.common.sm;

/**
 * Simplex.java
 *
 * Обобщенная реализация алгоритма симплекс-метода.
 * Для начальных данных использует интерфейс {@link ISMData}.
 *
 * Created on 11.11.2008 21:06:22
 * @author Дмитрий
 */
public class Simplex
{

  private ISMData smData;
  private int n;
  //базисная подматрица
  private double X[][];
  //обратная матрица
  private double X1[][];
  //целевая функция
  private double[] Res;
  //вектор двойственных переменных
  private double[] V;
  private double[] z;
  private double[][] X2;

  /**
   *
   * @param smData
   */
  public Simplex(ISMData smData)
  {
    this.smData=smData;
  }

  /**
   * Инициализация СМ
   */
  private void Initialization()
  {
    n=smData.getProblemSize();
    int i, j;
    X=new double[n][n];
    X1=new double[n][n];
    Res=new double[n];
    V=new double[n];

    //заполнение начального решения
    double t[][]=smData.getStartBaseMatrix();
    for(i=0;i<n;i++)
    {
      for(j=0;j<n;j++)
      {
        X[i][j]=t[i][j];
      }
    }
    //начальная обратная матрица и результат
    t=smData.getStartInverseMatrix();
    for(i=0;i<n;i++)
    {
      for(j=0;j<n;j++)
      {
        X1[i][j]=t[i][j];
      }
    }
    //начальный результат
    for(i=0;i<n;i++)
    {
      Res[i]=smData.getStartObjectiveFunction()[i];
      //smData.SMRes[n]+=smData.SMRes[i];
    }
    //вычисление вектора двойственных перменных
    //V[n]=1[n]*X1[n,n]   [V[M]=C[N']*a1[N',M]]
    for(i=0;i<n;i++)
    {
      V[i]=smData.getStartDualVariables()[i];
    }
  }

  /**
   * Проведение итераций СМ
   * @return
   * @throws trim.common.sm.SMException
   */
  private int Iterations() throws SMException
  {
    z=new double[n];
    X2=new double[n][n];
    SMColumn col;
    int icount=0;
    while((col=smData.getBestColumn(V)).getEstimation()>1e-11)//пока решение не отимально
    //while((epsilon=Maximize(r)+V[n]-1)>1e-11)
    {
      icount++;
      NextBase(col);//переход к новому базису
    }
    System.out.println(this.getClass().getSimpleName()+": SM iteration count="+icount);
    return 0;
  }

  /**
   * Преход к новому базису (итерация СМ)
   * @param col вводимый в базис столбец
   * @throws SMException
   */
  public void NextBase(SMColumn col) throws SMException//,int r[],double X2[][])//false-переход невозможен
  {
    int i, j, mx;
    double e;
    //double rrr;

    //новый план ищем в виде res[N]=res[N]-e*z[N]
    //z=X1[n',n]*r[n]   [Xt[N]=X[N]-t*z[N],t>=0]
    for(i=0;i<n;i++)
    {
      z[i]=0;
      for(j=0;j<n;j++)
      {
        z[i]=z[i]+X1[i][j]*col.getColumn()[j];
      }
    }    //if(mx==-1)
    {
      //поиск неотрицательного элемента
      mx=0;
      while((mx<n)&&(z[mx]<=1e-10))
      {
        mx++;
      }
      if(mx>n-1)
      {   //writeln(f,'Косяк');
        throw new SMException("Косяк");
      }
      for(i=mx+1;i<n;i++)
      {
        if((z[i]>1e-10)&&(Res[mx]/z[mx]>Res[i]/z[i]))
        {
          mx=i;
        }
      }
    }
    e=Res[mx]/z[mx];
    //уменьшение целевой функции
    for(i=0;i<n;i++)
    {
      Res[i]=Res[i]-e*z[i];
    }
    Res[mx]=e;

    //вводим новый столбец в базисную матрицу
    for(i=0;i<n;i++)
    {
      X[i][mx]=col.getColumn()[i];
    }

    //пересчет обратной матрицы
    //X1[n,n](новая)=D[n,n]*X1[n,n](старая) D[M,j]=E[M,j],j<>mx
    //D[n,mx]=(-z[1]/z[mx],..,1/z[mx],..,-z[n]/z[mx]
    for(i=0;i<n;i++)
    {
      for(j=0;j<n;j++)
      {
        if(i!=mx)
        {
          X2[i][j]=X1[i][j]-X1[mx][j]*z[i]/z[mx];
        }else
        {
          X2[mx][j]=X1[mx][j]/z[mx];
        }
      }
    }
    for(i=0;i<n;i++)
    {
      for(j=0;j<n;j++)
      {
        X1[i][j]=X2[i][j];
      }
    }

    //пересчет вектора двойственных переменных
    if(col.getEstimation()>0)
    {
      for(i=0;i<n;i++)
      {
        V[i]-=col.getEstimation()*X1[mx][i];
      }
    }
    /*for(i=0;i<n;i++)
    {
    V[i]=0;
    for(j=0;j<n;j++)
    V[i]=V[i]+X1[j][i];
    }*/
  }

  /**
   * Тестирование вводимого столбца(итерация СМ)
   * @param col вводимый в базис столбец
   * @return полученные коэффициенты целевой функци
   * @throws SMException
   */
  public double[] testNextBase(SMColumn col) throws SMException//,int r[],double X2[][])//false-переход невозможен
  {
    int i, j, mx;
    double e;
    //double rrr;

    //новый план ищем в виде res[N]=res[N]-e*z[N]
    //z=X1[n',n]*r[n]   [Xt[N]=X[N]-t*z[N],t>=0]
    double[] column=col.getColumn();
    for(i=0;i<n;i++)
    {
      z[i]=0;
      for(j=0;j<n;j++)
      {
        z[i]=z[i]+X1[i][j]*column[j];
      }
    }    //if(mx==-1)
    {
      //поиск неотрицательного элемента
      mx=0;
      while((mx<n)&&(z[mx]<=1e-10))
      {
        mx++;
      }
      if(mx>n-1)
      {   //writeln(f,'Косяк');
        throw new SMException("Косяк");
      }
      for(i=mx+1;i<n;i++)
      {
        if((z[i]>1e-10)&&(Res[mx]/z[mx]>Res[i]/z[i]))
        {
          mx=i;
        }
      }
    }
    e=Res[mx]/z[mx];
    //уменьшение целевой функции
    double[] resf=new double[n];
    for(i=0;i<n;i++)
    {
      resf[i]=Res[i]-e*z[i];
    }
    resf[mx]=e;
    return resf;
  }

  /**
   * Запуск итераций СМ
   */
  public void solve() throws SMException
  {
    Initialization();
    Iterations();
  }

  /**
   * Коэффициенты целевой функции
   * @return Коэффициенты целевой функции
   */
  public double[] getObjectiveFunction()
  {
    return Res;
  }

  /**
   * Значение целевой функции
   * @return Значение целевой функции
   */
  public double getObjectiveFunctionValue()
  {
    double r=0.0;
    for(int i=0;i<n;i++)
    {
      r+=Res[i];
    }
    return r;
  }

  /**
   * Вектор двойственных переменных
   * @return Начальный вектор двойственных переменных
   */
  public double[] getDualVariables()
  {
    return V;
  }

  /**
   * Возвращает базисную подматрицу
   * @return базисная подматрица
   */
  public double[][] getBaseMatrix()
  {
    return X;
  }

  public double[][] getReverseMatrix()
  {
    return X1;
  }

  /**
   * Возвращает исходные данные
   * @return исходные данные
   */
  public ISMData getSmData()
  {
    return smData;
  }
}
